Copyright © 1972 American Telephone and Telegraph Company 

The Bell System Technical Journal 

Vol. 51, No. 8, Octohcr, 1972 

Printed in U.S.A. 



Optical Power Flow in Multimode Fibers 

By D. GLOGE 

(Manuscript received May 8, 1972) 

Loss, coupling, and delay differences among the modes of multimode fibers 
influence their transmission characteristic in a complicated way. An 
approximation of the modes by a continuum leads to a comprehensive 
description of these interrelations. We relate the mode power distribution 
to the far-field output and calculate these distributions as functions of the 
fiber length and the input. We report measurements of the far-field dis- 
tributions at various lengths of a cladded low-loss multimode fiber. A 
comparison of theory and experiment yields a quantitative estimate of the 
mode coupling involved. We associate this coupling with random irreg- 
ularities of the fiber configuration and straightness, and construct a quanti- 
tative model of such irregularities. 

I. INTRODUCTION 

Some sources considered for use in optical communication systems 
have a spatially incoherent or multimode output and require overmoded 
fibers for efficient transmission. The fibers consist of a highly trans- 
parent core surrounded by a cladding of lower refractive index. Liquid 
core prototypes with losses as low as 20 dB/km have been built. 1 Solid 
mutimode fibers have slightly higher losses. 2 A recent study of their 
propagation and dispersion characteristics 3 showed a rather intricate 
behavior complicated by the fact that hundreds of modes could propa- 
gate simultaneously. These modes underwent a perpetual mixing process. 
The attenuation coefficient appeared to vary from mode to mode causing 
a relatively fast loss of the high-order modes. 4 An increase of delay with 
mode number (and fiber length) was observed as expected, but mixing 
and attenuation seemed to influence this relationship in a complex way. 

An exact knowledge of the processes involved is of considerable 
interest not only to understand the sources of loss in the fiber, but in 
order to determine the signal distortion in long fibers. It has been 
predicted 5 that under certain circumstances increased mixing reduces 
the signal distortion (ultimately forcing all energy to propagate at an 
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average velocity). But it remains to be determined what actual signal 
improvements can be gained in practice from this effect. Previous 
investigations of these problems, 5,6 although suited to show the concepts 
involved, were limited to model studies involving relatively few modes. 

In this work, we replace the modes by a continuum. This results in 
a relatively simple differential equation which describes the power 
distribution as a function of time, fiber length, and the continuous mode 
parameter. The differential equation can be solved rigorously for certain 
conditions which satisfactorily match experimental results. Explicit 
relations result which describe the propagation characteristics as a 
function of the modal coupling, attenuation, and delay coefficients. 
The coupling is then related to specific imperfections of the configuration 
or straightness of the fiber. 

This paper is primarily devoted to the time-independent solution of 
the problem. Signal distortion and, specifically, the (baseband) impulse 
response of long fibers can be derived from a slight modification of 
the above equations and this will be done in a subsequent paper. The 
concept underlying our results developed from experiments with solid- 
core fibers 3 but, in the meantime, measurements of long liquid-core 

fibers 7 have proven that these fibers follow the same concept. 

■ 

II. TRANSITION TO MODAL CONTINUUM 

For large mode numbers, the characteristic mode parameters change 
so little between neighboring modes that their discrete values can be 
replaced by one continuous variable. Consider the two-dimensional 
dielectric guide-a thin film, for example-sketched in Fig. 1. We assume 
that the relative index difference 



n 



(1) 



between core index n and cladding index n c is small compared to unity. 
In that case, the critical angle* for total internal reflection 



= aRS~^ 



(2) 



is small as well, and we can use small-angle approximations in the 
following relations. 

Within the high-index material, the field distribution of the rath 
mode is essentially sinusoidal (see Fig. 1) with transverse wave number 

* Defined here as the angle measured from the reflecting surface (see Fig. 1 ). 
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Fig. 1 — Sketch to illustrate the wave nature of the modes in the dielectric slab and 
in the fiber. 



U = 



inn 
2a 



(3) 



where 2a is the guide width. If A; = 2tt/X is the free-space wave number, 
the propagation direction of a mode (i.e., its representative plane wave) 
follows from 



nk 4an 

Because of Snell's law, this angle becomes 

" oui " 4a 



(4) 



(5) 



outside of the guide. In the far field (or the focal plane of a lens), the 
plane waves are concentrated about the directions +0 ou i and — out . 
The aperture of the guide determines the angular concentration of the 
two far-field "spots." If the guide width is 2a, the spot width is of the 
order of X/a. 

As we learn from (4), the propagation directions of neighboring 
modes differ by 



Ad = 



4an 



(6) 



and hence by X/4a outside the guide. The modes thus form a partly 
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overlapping sequence of spots in the far field, ordered according to 
mode number. Consequently the far-field distribution represents a 
direct image of the modal power distribution. 

The transition to the modal continuum uses a continuous angle 
instead of the discrete values (4). In this way, we arrive at a continuum 
of plane waves which, in the following, will be represented by rays. 
The power distribution P(0), in this continuum, is obtained by replacing 
ou t by in the (average) far-field power distribution. 

The cylindrical configuration lacks part of the conceptual clarity 
associated with the plane-wave representation, but a formal similarity 
permits us to arrive at an equivalent ray model which is satisfactory for 
almost all problems related to multimode fibers. We refer again to 
Fig. 1, considering now a cylindrical core of radius a imbedded in 
cladding material. The modal field distributions are given by Bessel 
functions. In the case of a small index difference, there are degenerate 
mode pairs (HE l + 1 , Q and EHi- Ull ) whose transverse mode number u is 
determined by the gth root of the Bessel function 8 

J l {ua) = 0. (7) 

Here I is the azimuthal order number. Fig. 2, which lists a few low-order 
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Fig. 2— The order numbers of degenerate fiber modes plotted versus ua and the 
effective group number m = 2au/r. 



POWER FLOW IN FIBERS 1771 

roots, has the purpose of indicating how much the exact roots deviate 
from the approximation 

»-s(« + 5). (8) 

which is to be used in the following. The group with q -\- 1/2 = 7, for 
example, which has ua = 22 according to (8), is marked by vertical 
lines in Fig. 2. 

Most problems of interest in multimode fibers (coupling, scattering 
loss, delay) require only the transverse wave number u for a satisfactory 
description of each mode. Furthermore, u can be related to a propagation 
angle and a far-field angle out — in exact formal agreement with (4) 
and (5) — through an Hankel transformation of the mode field at the 
fiber end. This transformation shows that a mode of azimuthal order I 
produces I far-field spots located on a circle which is defined by the 
angle out of (5). Figure 1 illustrates the situation if viewed as a meri- 
dional cross section through a cylindrical configuration. 

These facts suggest a description of the cylindrical modes by a single 
mode number 

m = 2q + I. (9) 

Equations (1) through (6) then obtain for the cylindrical guide as they 
do for the slab. The important difference is hidden in the fact that 
m of (9) comprises a group of modes with different q and I. As can be 
seen from Fig. 2, the number of possible combinations for a given m is 
the nearest integer below m/2. As mentioned earlier, every combination 
of q and I represents two (degenerate) modes. Consequently, each m 
describes a group of (approximately) m modes. In the far-field pattern, 
this mode group covers an annular area of "radius" out and approxi- 
mate "width" X/a. 

The transition to the continuum again converts to a continuous 
parameter. But is now considered as a radial variable which covers 
the solid angle ird] . We have a conceptual model which consists of a 
continuum of rays within the cone ird 2 c , whereby the modal power 
distribution P(0) is obtained by replacing out in the (average) far-field 
power distribution by of (4). 

To compute the total number of modes, we determine the highest 
possible group number m c = 4a?(0 c /X by inserting e into (4). If we 
consider also that each group has m modes; and each mode has two 
possible states of polarization, we have for the mode volume 

t 2 m = (*f )°2A. (10) 
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A comparison with the more accurate number (2iran/\) 2 A from Ref. 8 
gives an indication of the quality of the approximations used here. 

III. POWER FLOW EQUATION 

For the sake of simplicity, the following derivation is based on a 
model which seems to have limited validity at first glance. We assume 
that mode coupling takes place only between next neighbors. It will 
become apparent later that the error involved in this approximation is 
small if other modes couple also, but the coupling strength decreases 
sufficiently fast with the mode spacing. There is experimental evidence 3 ' 7 
that a mechanism of this kind is indeed present in real multimode fibers. 
Fig. 3, for example, shows a measurement performed with the solid-core 
multimode fiber mentioned previously. 2 ' 3 The core diameter was 55 /tm, 
the relative index difference A = 0.0046, and the nominal loss 33 dB/km. 
About 700 modes could propagate. By injecting a very narrow cone of 
light (through an index-matching cell), we excited about 150 of the 
low-order modes. This was measured by scanning the far-field of the 
output after 30 cm of fiber. Similar measurements with longer fibers 
revealed a slow but steady increase in the number of excited modes 
(the angular far-field width) with fiber length. This slow increase is 
considered as strong evidence of a power exchange which favors near 
neighbors and decreases rapidly with mode spacing. 

To simplify matters, let us again consider the two-dimensional case 
first. As long as the coupling mechanism is a statistical process, we can 
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Fig. 3— The increase of mode volume with fiber length measured in a multimode 
fiber which propagated 700 modes. 
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ignore the individual mode fields and obtain the power distribution 
directly from some form of power rate equations. 9 These consider the 
variation dP m in the power P m of the mih mode along a guide increment 
dz. In the time-invariant case, the variation dP m has two causes: (i) dis- 
sipation and loss (scattering) to the outside, which we comprise in a 
term — a m P m dz\ (ii) coupling to other modes. Our simplified model 
assumes coupling between neighboring modes only. Thus, if d m is the 
coupling coefficient between the modes of order m + 1 and m, we have 

^ = -« m P m + d m (P m + l - P m ) + Cl(P-l " Pm)- (11) 

dz 

The transition to the continuum requires power differences to be 
replaced by differentials. Especially, we set 

g=±i - P ~ = a. (12) 

m+1 - e m dd K J 

With 0,„ — m _i = A0 from (6), we can rewrite (11) in the form 

dPm D , A J, dP m , dP m -t \ , 1Q v 

— = -- m p m + w{d m -^ - 4^, ^rj • d3) 

The remaining difference requires an analogous transition 

, dP m , dP n -t Aa d (, dP m \ ,... 

After replacing the index m everywhere by a functional dependence of 0, 
we finally obtain the power flow equation 

f.-^ + ^i[««} (15) 

In the cylindrical case, the index m stands for a group of m modes. 
To obtain the power equation for the with mode group, we must there- 
fore sum (11) over all m members. The coefficients a m and d m depend 
only on m, and hence are the same for all group members. 10 However, 
the coupling to the lower group (m — 1) can occur only between m — 1 
members. Thus 

m ^t= = ~ma m P m + md m (P m+1 - P m ) + (m - l)d m _ 1 (P m _ 1 - P m ). 
dz 

(16) 
Using (12) and a transition analogous to (14), we obtain 
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With the help of (4) this finally leads to 

Because of the symmetry involved, we can expand a in the form 
a (0) = a + A0 2 + 

where a„ comprises loss common to all modes. A loss of this kind can 
later be accounted for by multiplying the final solution by a term 
exp ( — a z). For the moment, we ignore this part of the loss. Among 
the higher orders, the term Ad 2 is the most important one, because it 
essentially comprises the loss caused at the core-cladding interface. 
This is so because the power density at the interface increases quad- 
ratically with the transverse wave number u of a certain mode 8 and 
hence quadratically with 0. In the following, we retain only this im- 
portant term. 

The coupling coefficient d(9) requires the same expansion. Its zero- 
order term is essential and cannot be accounted for later on. Although 
no estimates exist on the magnitude of other terms, the following 
derivation merely retains this first term. Its physical significance will 
become clearer as we proceed. Thus with 

d(6) = d (19) 



or 



d = w d - = (LY d - • (2o) 



we can write (15) and (18) in the form 



^T = -Ad 2 P + D C§ for the slab, 

dz ad 



(21) 



and 



f =- A " 2p + 7^("^) ""• the fiber. (22) 



The form of the last terms in (21) and (22) identifies next-neighbor mode 
coupling as a diffusion process in the continuum. 

Solutions that are independent of z can be obtained from the sub- 
stitution 
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P = Qe~ yz , (23) 

where y denotes a (power) attenuation constant related to the "steady- 
state" solution Q. Equations (21) and (22) take the form 

D^H = {Ad 2 - y)Q for the slab, (24) 



30* 



and 



I Ye { e fe) = {A e " y)Q for the fiber ' (25) 



The first of these equations is satisfied by the Hermite-Gaussian, and 
the second by the Laguerre-Gaussian polynomials; both are well known 
from the theory of the open resonator. The attenuation parameters y 
associated with each of these solutions increase with the order of the 
polynomial. 

Both for the slab and the fiber, the solutions of least loss have the 
form 

exp (-6 2 /eZ) (26) 

with 

e„ = (4D/A)*. (27) 

The power loss associated with this distribution is 

Too = (AD) h for the slab, (28) 

and 

7„ = 2{AD) h for the fiber. (29) 

The distribution (26) constitutes an optimum balance between the loss 
in high-order modes and the steady outflow of power into those modes 
through coupling. It is assumed, of course, that the critical angle 8 e is 
so large compared to G M , that the steady-state distribution (26) is not 
significantly influenced by the boundary relations at = 6 C . If this 
is not the case, the solutions of (24) and (25) have to take these boundary 
relations into account. 

IV. BUILD-UP FROM GAUSSIAN INPUT 

Any z-dependent solution of (21) or (22) can of course be constructed 
from the infinite set of solutions of (24) and (25). But in the case of the 
fiber, there is a certain interest in special solutions which have an 
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arbitrary Gaussian input 

P ia = P exp [-07 62] (30) 

as initial condition. This is because it is convenient to study multimode 
fibers by using a Gaussian laser beam for excitation. A high-power lens 
converts this beam into the angular Gaussian distribution (30). By 
observing the change in P(6) with fiber length, and the loss as a function 
of various (Gaussian) input distributions, one obtains valuable infor- 
mation on the power flow in the fiber. 

Since both the input and the steady-state are Gaussian, it is reason- 
able to try the solution 

P - /(*) exp [-?/&*(?)]. (3D 

Although this approach is useful both for the two- and the three- 
dimensional configuration, we shall concentrate in the following on the 
fiber only. Introducing (31) into (22) yields the two differential equations 

dQ/dz = -^ 3 + 2D/6 (32) 

and 

df/dz = -iDf/Q 2 . (33) 

We can solve the first of these equations for 9 and obtain 

02 = e Jtanh 7eo (*+ 2o )| (34) 

(coth 7„(2 + z )) 

with the steady-state parameters e M and y„ from (27) and (29). The 
choice of tanh or coth and the coefficient z are determined by the 
initial conditions. Eq. (33) can be solved with the help of (34) and yields 

[cosh 7„(2 + Z )\ 

A Gaussian input stays indeed Gaussian, its width approaching monoto- 
nically that of the steady state. The transition function is the hyperbolic 
tangent if the input width is smaller than the steady-state width and the 
hyperbolic cotangent in the opposite case. For the initial conditions (30), 
the solutions (34) and (35) can be written in the form 

6 2 (e) = G& g + gtanh 7 ,s 
0«, + Of tanh 7 M 2 
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and 

P Q 2 
I(Z) = elsmhyjs+'elcoOiyjB' (37) 

To obtain the total power in the guide, P(6) must be integrated over 
all angles up to 8 C . If we assume, as previously, that P(6) is sufficiently 
small at the critical angle C and beyond, we can extend the integration 
to infinity. With (31), the total power is 

2t f P(8)ddd = icfQ 2 . (38) 

Jo 

The power loss per unit length is consequently 

t»--j£i<W. (39) 

By using the differentials dQ/dz and dj/dz from (32) and (33), we obtain 

y(z) = AQ 2 (z). (40) 

With (27) and (29) this can also be written in the form 

7(2) e 2 ( 2 ) 



e 



(41) 



The ratio (41) is plotted in Fig. 4 versus the fiber length for some specific 
input conditions. The plot illustrates the loss, the solid angle covered 
by the fiber output and, since this is proportional to the mode volume, 
also the number of modes propagating in the fiber. 

If the measured width 6(2) is small compared to the steady-state 
width 6,,, , we can approximate (36) by 

e 2 = e^oo* + e 2 , . (42) 

In this case, because of (27) and (29), 2 (z) is a straight line with the 
slope 

e^7- = 4D. (43) 

Thus measuring 6 2 (z) under these conditions yields directly the coupling 
parameter D. The approximate linear increase of the data in Fig. 3 is 
an indication of the validity of (42). A straight-line approximation of 
the measured data yields 

D = 7 10- 5 rad 2 /m. (44) 

This value represents a first approximation for the zero-order coupling 
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Fig. 4 — Calculated increase of mode volume and loss with fiber length. The data 
are related to the steady-state values. 



coefficient d which, because of (20), becomes d = 16 m -1 at the measur- 
ing conditions (a = 50 /im, X = 0.63 /im). The range of the measured 
data is not sufficient to draw any conclusions on the size of higher-order 
coefficients in d(6). 

V. TWO SOURCES OF COUPLING 

Marcuse 10 has studied a dielectric slab guide with slightly distorted 
interfaces. He finds that two modes are coupled if the surface imper- 
fections comprise a component of wavelength A that coincides with 
the "beat wavelength" between the two modes. The beat wavelength 
is the distance in which the phase difference between two modes increases 
to 2ir. It can be calculated from the wave number 



(fcV - w 2 )* « kn - - j- 



(45) 



of the mth mode. With u from (3), we obtain for two neighboring modes 
of order m and m + 1, 
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2ir 4a fAa , 

A " KT^T, = T (46) 

This is also the distance that a ray of angle requires to complete a 
zigzag period in a slab of width 2a. 

To describe the distorted slab walls, let us define a "power" spectrum 
0(1 /A) of the (random) deviations 8(z) from ideally straight interfaces. 
If the irregularities of both walls are uncorrelated, the coupling coeffi- 
cient d(d) between adjacent modes is 10 

-lews- 

Using this relation and (20), we can calculate the spectrum associated 
with the measured coupling parameter D of (44): The result 

4) - §£ » 

suggests a decrease of <f> with the fourth power of the (spatial) frequency 
1/A. 

This result permits us to estimate the coupling among nonadjacent 
modes. As is evident from (45) and (46), modes which differ by a 
small number r have the beat wavelength 4ra/0. Because of (48), 
coupling among such modes decreases with the fourth power of their 
order difference. It is this rapid decrease which permitted us to neglect 
all coupling except that between next-neighbors in (11). The error 
involved in this approximation can be estimated for the case that the 
power distribution P(6) changes slowly within rAd. In that case, 

The transition from (11) to (15) then allows us to lump the coupling 
of all modes from m to m + r in the coefficient d m = d(d) which assumes 
the form 

dtp = ^d^l.0Sd. (50) 

This suggests that the relative error in our approximation is as small 
as 8 percent. 

Random bends in the guide are another source of coupling. This 
problem has been studied by tracing rays through a randomly curved 
slab guide. 11 The result relates the statistics of the ray angle d to the 
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"power spectrum" C(l/A) of the curvature components. Reference 11 
demonstrates that components of wavelength A predominantly influence 
rays with the same zigzag wavelength. For a short guide length z, the 
probability distribution of the ray angle 6 is found to have a variance 
which increases as 

where o- 2 is the variance at the input. 11 These results presuppose Gaussian 
statistics, for which the probability to find a ray at 6 has the form 
exp (-0 2 /2A 

Let us compare this distribution to the distribution 

P(0) = exp[-0 2 /e 2 ] 

of (31). For negligible mode attenuation, the width of this Gaussian 
is given by the simple relation (42). Like the variance of (51), it grows 
linearly with length. A comparison of the growth factors involved must 
take the factor 2 into account which enters because of the definition of 
the variance. This leads to the relation 

To compare this result with <t> of (48), let us consider correlated 
deviations of the form 5 sin (2tz/ A) at both walls. Twofold differentia- 
tion with respect to z transforms this into a curvature component of 
the form (27r/A) 2 5 sin (27rz/A). Accordingly, we can relate the curvature 
spectrum to a spectrum of (correlated) irregularities 

*.Q) = (A/27r) 4 C. (53) 

This result transforms (52) into 

The factor 2 which distinguishes this result from (48) results from the 
correlation of the wall deviations assumed in (53) contrary to (48) ; 10 to 
represent a curved slab, the deviations must be equal and in phase. 
Irregularities of this kind couple only modes which differ by an odd 
order number r = 1,3, • • • . Coupling across even numbers results from 
irregularities in anti-phase or — in the case of the fiber — from irregularities 
of at least twofold cross-sectional symmetry. Specific parts of the spec- 
trum 0(1/ A) are likely to be dominated by certain kinds of irregularities. 
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The region of interest is determined by A = 4a/6 and was in our case 
between A = 1 and 10 mm. Since these lengths are significantly larger 
than the core diameter (150 nm), we believe that bends were the domi- 
nant source of coupling. The following example therefore uses the 
relation (52) for a quantitative estimate of the irregularities involved. 
Although (52) applies specifically to the slab model, we shall combine 
it with the fiber data (44), confident that this will illustrate at least 
the orders of magnitude involved. 

To obtain a more tangible description of the random curvature, we 
assume it to be composed of randomly distributed singular deviations 
of the kind illustrated in Fig. 5a. We model these deviations by single- 
period sinewaves of the form 5 sin (2irz/A). Essentially, only those with 
a width larger than A/2 contribute to the curvature spectrum at A. 
If there are rj of those per unit length, the curvature spectrum has 
approximately the value 11 

C = St 4 V 5 2 /A 2 (55) 

in the vicinity of A. Because of (52) 

18 = 3^- (56) 

Fig. 5b evaluates this relation for the case of the fiber measured. Plotted 
is the density y versus the amplitude 8 for ray periods (and angles) of 
interest. For example, an average of 1000 singular irregularities per meter 
would account for the coupling measured, if their magnitudes obeyed 
the relation 8 = 15- 10 _0 A. That would mean that the magnitude of the 
irregularities increases linearly with their length reaching a value of 
15 nm at A = 1 mm. 

VI. CONCLUSIONS 

A comprehensive description of a multimode fiber by one differential 
equation is possible, if the modes are approximated by a continuum. 
Under certain realistic conditions, this equation has a rigorous solution. 
We calculate here the far-field output distribution as a function of fiber 
length and compare this to experiments performed with a low-loss 
multimode fiber. We find neighboring modes to be coupled by an average 
1.6 percent per mm. Coupling among other modes seems to be at least 
an order of magnitude less. 

Among the possible sources of coupling taken into consideration, we 
believe random bends to be the most likely. In this case, the curvature 
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Fig. 5 — Deviations from straightness as indicated in (a) account for the coupling 
measured, if their density jj, and magnitude 5 are related to their length A as plotted 
in (b). 



spectrum has a value of 2 • 10" 



in the region of spatial ferquencies 



between 0.1 and 1 mm -1 . Another equivalent description of this result 
is by small hump-shaped deviations from straightness — on the average 
1000 of them per meter— the magnitude of which increases propor- 
tionally to their length and is about 15 nm for a length of 1 mm. 

The theory derived here can be modified to include the velocity dif- 
ferences among the modes and, in this way, to describe the impulse 
response in the presence of coupling and loss. This will be the subject 
of another paper. 
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